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Abstract. It is shown how the Beneath-and-Beyond algorithm can be used to yield 
another proof of the equivalence of V- and H-representations of convex polytopes. 
In this sense this paper serves as the sketch of an introduction to polytope theory 
with a focus on algorithmic aspects. Moreover, computational results are presented 
to compare Beneath-and-Beyond to other convex hull implementations. 



1 Introduction 

One of the most prominent algorithmic problems in computational geometry 
is the convex hull problem: Given a finite set of points S* S IR*^, compute 
the facets (for instance, in terms of their normal vectors) of the convex hull 
conv S. To the non-expert it may come as a surprise that fundamental ques- 
tions concerning the complexity of this problem are still unsettled. Numerous 
methods have been invented through the last three decades. However, the per- 
formance of each known algorithm heavily depends on specific combinatorial 
properties of the input polytope conv S. Our current knowledge can largely 
be summarized by three statements: For many algorithms there is a class of 
polytopes for which the given algorithm performs well. For all known algo- 
rithms there is a class of polytopes for which the given algorithm performs 
badly. There are classes of polytopes for which all known algorithms perform 
badly. 

For a comprehensive survey on convex hull algorithms we refer the reader 
to the paper [3] by Avis, Bremner, and Seidel. 

The Beneath-and-Beyond algorithm is among the most natural methods 
for convex hull computation. It is treated thoroughly in many text books, 
see, for instance, Griinbaum [19, 5.2] and Edelsbrunner [9, Section 8.4]. The 
purpose of this paper is twofold: Firstly, we want to sketch how this partic- 
ular algorithm can be used to give yet another proof of the "Main theorem 
for polytopes", see Ziegler [30, Theorem 1.1], which says that the convex 
hull of finitely many points is exactly the same as a bounded intersection of 
finitely many affine halfspaces. Previously known proofs are based on Fourier- 
Motzkin elimination, as in Grotschel [18, Kapitel 3] and Ziegler [30], on the 
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simplex method for linear optimization, see Schrijver [28, Chapter 7], or non- 
algorithmic, as in Griinbaum [19, 3.1.1] or Matousek [25, Theorem 5.2.2]. Our 
strategy for a proof of the "Main theorem for polytopes" via the Beneath- 
and-Beyond algorithm could be phrased as: Prove everything directly for 
simplices and then inductively use triangulations to extend the results to 
arbitrary polytopes. Secondly, we give a brief description of the implementa- 
tion of the Beneath-and-Beyond algorithm in polymake [14-16]. The paper 
closes with a survey of computational results of Beneath-and-Bcyond in com- 
parison with Fourier-Motzkin elimination (or, dually, the double description 
method), implemented by Fukuda [12], and reverse search, implemented by 
Avis [2]. This complements the computational results in [3]. 

I am indebted to Thilo Schroder who helped to obtain the timing data in 
Section 5. Thanks to Volker Kaibel, Marc E. Pfetsch, and Giinter M. Ziegler 
for their comments and many helpful discussions. And, finally, thanks to Ew- 
genij Gawrilow for co-authoring polymake and his unfailing technical advice. 
The polytope images were produced with JavaView [27] and polymake. 

2 Definitions, an Algorithm, and a Classical Theorem 

A subset S of the Euclidean space is convex if for any two points x,y € S 
the line segment [x,y] = { Ax + (1 — A)y ] < A < 1} between x and y is 
contained in S. We define the convex hull of S, which is denoted by conv(S'), 
as the smallest convex set containing S. It is easy to see that, equivalently. 



For the purpose of this paper the key objects are (convex) polytopes, that 
is, the convex hulls of finite point sets in R'^. 

An affine hyperplanc H defines two (closed) affinc halfspaccs _ff + and H~ 
with H+UH- = R'^ and H+ D H' = H. Let 5 C R'^ be any set. A hyper- 
plane H with iJ+ n 5 ^ and H~ n S ^ is said to separate S. A hyper- 
plane H supports the set S if H intersects S non-trivially and if, moreover, 
either S C or 5 C H~. We will always assume that S C H~^, that is, if 
H supports S we consider H to be positively oriented toward S. 

The dimension of any set 5 C R'' is defined as the dimension of its affine 
span aff(S'); it is denoted by dim 5. 

Throughout the following let X c R'^ he finite. A proper face of the 
polytope P = conv(X) is the intersection of P with a supporting hyperplane 
which does not contain P. Note that this last condition is superfiuous in the 
case where dimP = d, that is, aff(P) = R'^. Occasionally, the empty set and 
the set P are also considered as (non-proper) faces of P. Some types of faces 
of a polytope will play particularly important roles in our investigations: the 
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faces of dimensions 0, 1, d — 2, and d — 1, are called vertices, edges, ridges, 
and facets, respectively. 

A point X G S' C R'' is a relatively interior point of S if there is a small 
dim 5-dimensional ball around x which is contained in S. The boundary of S 
is the complement of the relative interior of S in S; it is denoted by bnd(S'). 

For each facet of a polytope P one can choose a facet defining afHne 
halfspace whose boundary hyperplane supports P in F. A complete set 
of facet defining halfspaces contains exactly one such halfspace for each facet. 
If dim P = d then this choice is unique. 

In the following we abbreviate "rf-dimensional polytope" and "fc-dimen- 
sional face" by "d-polytope" and "fc-face" , respectively. The number of fc-faces 
of a polytope P is denoted as fk{P)- The sequence (/o(P), . . . , fd-i{P)) is 
called the /-vector of the d-polytope P. 

Consider n+1 affinely independent points xi, . . . , Xn+i € IR*^. Their convex 
hull is an n-simplex. A direct computation yields the following result. 

Proposition 1 The k-faces of any n-simplex An ~ conv(a;i, . . . ,Xn+i) are 
precisely the convex hulls of all k + 1 element subsets of {xi, . . . ,Xn+i}- In 
particular, the proper k-faces are k-simplices and fk{An) = {2+1) ■ Moreover, 

each n-sim,plex in R" is the intersection of its n+\ facet defining halfspaces, 
and the boundary of a simplex is the union of its facets. 

Observe that 



is an (n — 1) -dimensional simplex. Hence the Equation (1) implies that each 
polytope is the linear projection of a high-dimensional simplex. 

A (geometric) simplicial complex is a finite collection T of simplices in IR'' 
with the following two properties: 

a. Each proper face of a simplex in T is also contained in T. 

b. The intersection of any two simplices A, A' G T is a (possibly empty) 

face of both, A and A' . 

A simplicial complex T is a triangulation of a set 5* C R'' if the union of all 
simplices in T is 5. The fc-dimensional elements of a simplicial complex T 
again are called k-faces of T. A trivial example: Any simplex (together with 
its collection of faces) is a triangulation of itself. 

Throughout the following we will assume that simplicial complexes (and 
triangulations) are pure, that is, all (with respect to inclusion) maximal faces 
have the same dimension. 

Usually, it is more convenient to assume that a given polytope P = 
conv(Ar) affinely spans its ambient space R"^. This is justified — also in an 
algorithmic setting — due to the following reasoning. By performing Gaus- 
sian elimination we can determine the dimension dimP = dimX, and we can 
even select an affinc basis of the span of X . Moreover, for instance by omit- 
ting redundant coordinates, we can project X afhnely isomorphic to a linear 
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subspace L of R'' with dimL = dimX. A triangulation of the projection of P 
can directly be hfted back to P C R''. Similarly for facet defining halfspaces. 

Lemma 2 Let T he any triangulation of a d-polytope P c R'^. Then the 
following holds. 

a. For each face F of P the set T{F) = {A €:T \ A ^ F} is a triangulation 
ofF. 

b. For each {d — l)-face AofT contained in the boundary of P there is a 
unique facet of P which contains A. 

Proof. Let H he a hyperplanc which supports P such that the intersection 
HnP is the face F. We have to show that for each point x G F there is a face 
A € T with X € A and A C F. For each face of T the intersection with H 
is again a face of T. But, T covers P, that is, some face of T contains x. Its 
intersection with H is the desired face A. This proves the first statement. 

Now let A gT be a {d — l)-face in bnd(P). Since T is pure, there exists 
a d-face A' e T with the property that Zi is a facet of the d-simplex A'. 
Therefore, there is a hyperplane H with if+ D A' and H Ci A' = A. Choose 
a point y in the relative interior of A. Suppose that H separates P, that is, 
there is a point x S H~ n P. Observe that A" = conv(Z\, x) is a rf-simplex, 
which is contained in P since P is convex. Now A' U A" is a d-dimensional 
ball which contains y in its interior. This contradicts y G bnd(P), and thus 
H defines a facet of P. 

A direct consequence of the preceding lemma is the correctness of the 
Algorithm A below, which computes the complete set of facet defining half- 
spaces of a polytope from a given triangulation. 



Algorithm A: Extracting the facets of a polytope from a triangulation. 

Input : triangulation T of P = conv(X) 

Output: complete set of facet defining halfspaces of P 

foreach (d — l)-face A gT do 

if a,S{A) does not separate X then 
L J^U{aff(Z\) + } 

return 



Now we are ready to prove the main result of this section: Each polytope 
is a bounded intersection of finitely many affine halfspaces. Observe that the 
actual statement of the theorem is, in fact, much stronger. This is necessary 
to allow for an easy inductive proof. 
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Theorem 3 Let P = conv(X) be a polytope and T a complete set of facet 
defining affine halfspaces. Then the following holds. 

a. There is a triangulation T of P such that the vertices of T are precisely 
the points in X . 

b. The polytope P is the intersection of aff (P) with the intersection of all 
halfspaces in T . 

c. The boundary of P is the union of its facets. 

Proof. We give a constructive proof. More precisely, starting from the finite 
set X, we construct T and T with the desired properties. As pointed out 
above we can assume that P is fuU-dimensional. 

For the rest of the proof we fix an arbitrary ordering of the set X = 
{xi , . . . ,Xn} such that the first d+1 points xi, . . . , Xd+i are affinely indepen- 
dent. This ordering gives us a sequence of (i-polytopcs Pk = conv{a;i, . . . , Xk} 
for k e ... , n}. We have Pk+i = co-nv{Pk,Xk+i) and we proceed by 

induction on k. 

Since Xi, . . . ,Xd+i are affinely independent, their convex hull Pd+i is a 
simplex. The vertices xi, . . . ,Xd+i of Pd+i clearly are the vertices of the 
trivial triangulation of Pd+i by itself. The facet defining halfspaces J-d+i of 
Pd+i are given by Proposition 1. And we have that pl^^+i = Pd+i and Pd+i 
is the union of its facets. 

For the inductive step suppose that we have a triangulation 7^ and the 
facet defining halfspaces J^k of Pfe. We construct a triangulation Tk+i by 
using Algorithm B below, where we let P = Pk and a; = a;fe+i. It is clear that 
this yields a simplicial complex, but we have to prove that lj7fe+i is Pk+i. 
Without loss of generality Xk+i ^ Pk, otherwise Pk+i = Pk and the algorithm 
does not change the triangulation. 

The induction hypothesis provides us with a facet F of Pk which is vio- 
lated by Xk+i ■ By Lemma 2a we have a. {d— l)-face AoiTk which is contained 
in F, and conv(Z\, a;fc+i) G Further, consider some point x G Pk. By 

compactness, the line segment [x, Xk\ meets the boundary of Pk in a point. 
Due to compactness the line segment [a;,a;fe+i] meets the boundary in a (not 
necessarily unique) point y. Again by induction, bnd(Pfe) is the union of 
its facets, so y is contained in some facet of Pk and therefore also in some 
{d — l)-face of 7^. We obtain 



and hence Tk+i is a triangulation of Pk+i- 

Now the Algorithm A can be used to determine the complete set J^k+i 
of facet defining halfspaces of Pk+i- It is clear that Pk+i C f]J^k+i- For the 
reverse inclusion suppose that z is a point in f^J-'k+i \ Pfc+i. By repeating 
the same argument as before we obtain a (d— l)-face A e Tk+i in bnd(Pfc+i) 
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such that a.S{A) is a hyperplane which intersects Pk+i in a facet and which 
separates z from Pfc+i- This yields the desired contradiction. 

Making use of both parts of Lemma 2 allows to conclude that bnd(P;j+i) 
is the union of the facets of P^+i. The theorem is proved. 



Algorithm B: Extending a triangulation. 

Input : triangulation T of rf-polytope P, facet defining halfspaces 

of P, point X 
Output: triangulation of P' = conv(P, x) 

foreach F+ G do 
if X <^ F+ then 

foreach {d - l)-/ace Z\ G T with A C bnd(P+) do 
|_ T' <- T' U {faces of conv(zi, x)} 

return T' 



As a direct application we obtain Caratheodory's Theorem: 

Corollary 4 Let X G R''' be a (full- dimensional) finite set of points and p G 
conv(X). Then there is an affinely independent subset X' = {xq, . . . C X 
with p G conv(X'). 

Proof. By Theorem 3 the polytope conv(X) admits a triangulation T with 
vertex set X. For X' choose the set of vertices of any d-simplex in T which 
contains p. 

This is where we end our algorithmically inspired introduction to polytope 
theory. Of course, this is not the end of the story. The next step, which then 
paves the way to the rest of the theory, would be to prove the converse of 
Theorem 3: Each bounded intersection of finitely many afiine halfspaces is a 
polytope. One classical way of proving this is via the Separation Theorem, 
see Matousek [25, Theorem 1.2.4], and duality, see [25, 5.1]. Then Theorem 3 
itself can be applied to prove its converse. It is interesting that it is also 
possible to reverse the order in which one proves these results. 

3 Sizes of Triangulations and Algorithm Complexity 

One of the major open questions in computational geometry is whether there 
is a polynomial total time convex hull algorithm, that is, an algorithm whose 
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running time is bounded by a polynomial in the number of vertices and facets. 
It is unreasonable to hope for an algorithm whose running time depends 
polynomially on the input size only, since there are families of polytopes 
whose numbers of facets grow exponentially with the number of vertices. For 
instance, take the rf-dimensional cross polytopes 

conv{±ei, . . . , ±ed} 

with 2d vertices and 2*^ facets. 

In the following we relate the complexity of Beneath-and-Beyond to the 
size of the triangulations produced. Then we review some explicit construc- 
tions of polytopes with large triangulations. 

3.1 Complexity Analysis 

Wc want to examine how the Beneath-and-Beyond method fits into the pic- 
ture: What is its complexity? Rather than delving into technicalities we want 
to exhibit the geometric core of this question. The very coarse and schematic 
description of the algorithm in steps A and B overestimates the costs of the 
Beneath-and-Beyond method. For hints to a more practical approach see the 
next section. 

We begin with Algorithm A, which extracts the facets of a rf-polytope P 
from a triangulation T. We call the number of d-faces of T the size of T. 
li t = size(T) then T has at most {d + l)t faces of dimension d — 1. A 
facet normal vector can be computed from d affinely independent points con- 
tained in the facet by Gaussian elimination, which requires 0{d^) steps. The 
overall complexity of Algorithm A is then bounded by 0{d^nt), where n is 
the number of vertices of P. Now consider Algorithm B: The input is a d- 
polytope Pk with ruk facets and Uk vertices, a triangulation 7^ with tk faces 
of dimension d, and one extra point Xk+i- The desired output is a trian- 
gulation of Pfe+i = conv(Pfe, Xfc+i), which can be computed by evalutaing 
0{dmktk) scalar products since each d-simplex of 7^ contains exactly d -\- \ 
simplices of dimension d—1; this gives a total of 0{d^mktk) arithmetic oper- 
ations. Summing up we obtain an upper bound for the complexity of a single 
Beneath-and-Beyond step. 

Lemma 5 The vertices, the facets and the triangulation T^+i of Pu+i can 
he computed from the triangulation Tk of Pk in 0{d'^ m.a,x{mk,nk)tk) steps. 

We have tk < ifc+i and ruk < (d + l)ife. Setting m = m,n and t = tn we 
can sum up for all A; G {d -|- 1, . . . , n} to obtain the following result. 

Proposition 6 The overall complexity of the Beneath-and-Beyond algorithm 
is bounded by 0{d^nt^). 
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Of course, an implementation of the algorithm as sketched in Algorithms A 
and B is far from optimal, see Section 4 below for a few more details. More- 
over, the analysis given is very coarse and could be sharpened easily. However, 
for the purpose intended it is good enough: Since the size of the final trian- 
gulation clearly is a lower bound, the Beneath-and-Bcyond algorithm runs in 
polynomial total time if and only if the size of the triangulation constructed 
is bounded by a polynomial in m and n. 

And this could be the end of the story, because it is known that products 
of simplices form a family of polytopes where the size of any triangulation 
is super-polynomial in the number of vertices and facets, see Haiman [20] 
or Avis, Bremner, and Seidel [3, Lemma 3]. We conclude that Beneath-and- 
Beyond has a worst case super-exponential total running-time. 

It is a consequence of the Upper Bound Theorem, see Ziegler [30, Sec- 
tion 8.4], that the parameter t is bounded by 0(n'-'^/^-l ). This bound is ac- 
tually attained, for instance, by the cyclic polytopes, see [30, Example 0.6], 
which arise as convex hulls of finitely many points on the moment curve 

t^{t,t\t\...,t''). 

If the input points are sorted it is possible to avoid looking at all the facets 
of Pfc in Algorithm B. In this way Algorithm B can be replaced by a method 
which takes time which is proportional to the number of facets of Pk which 
are not facets of Pk+i- Amortized analysis then shows that the size t of the 
final triangulation 7^ only enters linearly into the cost function. In particular, 
for fixed dimension d, one obtains an 0(nlogn -|- ri'-(''+^)/^-l ) algorithm, see 
Edelsbrunner [9, Section 8.4.5]. 

There is a fairly general result due to Bremner [6] who proves that each 
incremental convex hull algorithm has a worst-case super-polynomial total 
running-time, where incremental means that the algorithm has to compute 
the convex hulls of all intermediate polj^opes Pk- In particular, this also 
proves that Beneath-and-Beyond is not a polynomial total time algorithm. 

3.2 Polytopes with Large Triangulations 

In spite of the fact that there is no polynomial total time convex hull algo- 
rithm known, some of the known algorithms have a polynomially bounded 
running-time on special classes of polytopes. In particular, the reverse search 
method by Avis and Pukuda [4] runs in 0{dmn) time on simplicial poly- 
topes. Each simplicial polytopc has a small triangulation: Choose a vortex v 
and cone (with apex v) over all the facets not passing through v. Such a 
triangulation is extremely small, it is clearly of size 0{m). Since the size of a 
triangulation is the decisive factor for Beneath-and-Beyond's time complex- 
ity, this raises the question whether Beneath-and-Beyond could at least have 
a polynomial total running-time on simplicial polytopes. We note that the 
polytopes in Bremner's construction [6] arc products of simplicial polj^topes, 
that is, they are neither simplicial nor simple. 
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A placing triangulation of a polytope P is a triangulation produced by the 
Beneath-and-Beyond algorithm for some ordering of the vertices. For a given 
ordering, the size of the corresponding placing triangulation is related to the 
question of how many facets of the intermediate polytope Pk are violated by 
the next vertex Vk+i ■ A placing triangulation with respect to a given vertex 
order is the same as a pushing triangulation with respect to the reverse order, 
see Lee [24]. In particular, placing triangulations are lexicographic and thus 
regular. 

In the dual setting, that is with the roles of vertices and facets inter- 
changed. Avis, Brcmncr, and Scidcl [3] introduced the concept of dwarfing: 
Informally speaking, a separating hypcrplanc H is called dwarfing for a poly- 
tope P if all facets of P arc also facets of P fl , and if very many vertices 
of P are not vertices of P Pi . Equivalently, the star of the vertex vh dual 
to H in any placing triangulation of the polar polytope (Pni?+)*, where Vh 
comes last, is large. 

The regular d-dimensional cube Cd = [0, 1]'' has 2'' vertices and 2d facets. 
Consider the affine halfspace = {a; G IR'' | — ^/^}- The boundary 
hyperplane Hd = dH^ separates the cube Cd- Now, as in [3, Theorem 4], call 
the simple d-polytope 

Cd^CdH H+ 

the dwarfed d-cube. It has 2d+l facets and d"^ + 1 vertices: The origin and 
the d unit vectors ei , . . . , are the only vertices of Cd which are contained 
in H^; there is one new vertex + ^ej for each i, j G {1, . . . , d} with i ^ j. 
For illustrations see Figures 1 and 2. By Barnette's Lower Bound Theorem, 
see Br0nsted [7, §19], rf^ + 1 is the minimal number of vertices for a simple d- 
polytope with 2d+ 1 facets. It follows that the dwarfed cubes can be obtained 
from a simplex by repeated truncation of faces. 

Since the cube and the dwarfed cube have d+1 vertices in common, there 
are 2'* — d — 1 vertices of Cd which are not vertices of Cd- We conclude the 
following result. 

Proposition 7 Let vi, . . . , V2d+i be an ordering of the vertices of the polar 

Cd of the dwarfed d-cube Cd such that the last vertex v' = V2d+i corresponds to 
the dwarfing facet Hd ofcd- Then the number of d-simplices which contain v' 
in the induced placing triangulation of c*^ equals 2*^ — ^—1. 

The dwarfed cubes form a family of polytopes that are "bad" as input for 
the Beneath-and-Beyond algorithm. Since there is only one dwarfed cube per 
dimension, this docs not tell anything about the situation in fixed dimension. 
However, the same idea can be applied to another (bi-parametric) family 
of polytopes, whose number of vertices and facets is unbounded even if the 
dimension is fixed. We sketch the construction from [3], and we omit the 
proofs. 
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Fig. 1. Dwarfed 3-cube C3 with dwarfing facet marked. 



For d = 26 > A and s > 3 let Gd{s) be the d-dimensional polytope defined 
by the foUowing hst of 6s hnear inequalities, all of which define facets: 



SXk 

{2i + l)xk 
(2s - 3)xk 



yk>0 
-yk>0 

-Vk < i2i + l){s + i) 
-Vk < 2s(2s-3), 



(2) 
(3) 
(4) 
(5) 



where i G {0, . . . , s — 4}, k G {1,...,6}, and a vector in R'* is written 
as (xi, j/i, 2:2, ?/2, • ■ • tXSiVs)- The polytope Gd{s) is the product of 6 copies 
of the s-gon obtained by fixing k in the list of inequalities 2 to 5. This 
product of s-gons has vertices. There is a dwarfing halfspace H^^ = 
{x £ R'^ \ ^ < 2s — 1}, and the dwarfed product of s-gons 



gd{s) = Gd(s) n 



is simple. It has (5s + 1 facets but only {d — l){6s + 1 — c?) + 2 vertices, which, 
again, meets Barnette's lower bound; see [3, Theorem 6]. 

Since Gd{s) and gd{s) can have at most {d — l)((5s + 1 — c?) + 1 vertices in 
common, this yields a similar result as for the dwarfed cubes. 
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Fig. 2. Explosion of the placing triangulation of the polar polytope C3 induced 
by the ordering indicated as vertex labels. Only the four 3-simplices containing 
the final vertex (numbered 7), which corresponds to the dwarfing facet of C3, are 
displayed in solid. 



Proposition 8 Let Vi, . . . , f (ci~i)(5s+i-d)+2 be an ordering of the vertices of 
the polar gd{s)* of the dwarfed product of polygons gdis) such that the last 
vertex v' — W((i-i)((5s+i-d)+2 corresponds to the dwarfing facet Hd^s of gd{s). 
Then the number of d-simplices which contain v' in the induced placing tri- 
angulation of gd{s)* is at least — {d — l)(^s + 1 — d) — 1. 

This already shows that the worst case running time of Beneath-and- 
Beyond is not polynomiaUy bounded in combined the size of the input and 
the output total time algorithm. Additionally, one can show that Gd{s) and 
gdis) share exactly 6{s — 2) -(- 1 vertices, so the precise number of d-simplices 
containing v' is — 6s + 26 — 1. 

The closer analysis in [3, Theorem 12] reveals that essentially from the pre- 
ceding Proposition 8 it follows that even a typical placing triangulation (that 
is, with respect to a random ordering) of gd{s)* grows super-polynomially 
with the number of vertices and facets. 

Theorem 9 The polar dwarfed products of polygons gd{s)* «re simplicial d- 
polytopes with (d — l){ds/2 + 1 — c?) + 2 vertices and ds/2 + 1 facets, but with 
an expected size of a placing triangulation of order Qi^s'''^'^ / {d + 1)). 
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It seems to be open whether there is a class of simpHcial polytopes such 
that each placing triangulation is large compared to their number of (vertices 
and) facets. 

4 On the Implementation 

In the preceding section we related the performance of the Bcneath-and- 
Beyond algorithm to the size of certain triangulations. While this captures the 
main ideas of the method, an implementation which is feasible for practical 
problems is slightly more involved. We sketch the well-known key points 
below; for a more thorough discussion see Edclsbrunncr [9, 8.4.5]. 

Rather than extracting the facets of the intermediate polytope Pk from the 
triangulation Tk from scratch (as Algorithm A suggests) it is more natural to 
store the facets and their neighborhood structure in terms of the dual graph 
of Pk- Adding the next vertex Vk+i then requires: (i) to find one violated 
facet, (ii) to perform a breadth- first-search in the dual graph of Pk to find 
all the other violated facets, (iii) to extend the triangulation Tk of Pk to the 
triangulation T^+i of Pk+i by coning over the induced triangulation of the 
violated facets, (iv) to determine the new facets, that is, the facets of Pk+i 
which are not facets of Pk by examining Tk+i- 

Edelsbrunner [9, 8.4.5] advocates to sort the input points lexicographically 
in order to eliminate the time for finding the first violated facet, which then 
allows to find all the violated facets in linear time (in fixed dimension). Our 
implementation does not rely on ordered input, but it is programmed in a 
way such that Edelsbrunner's analysis applies, if the input happens to be 
ordered. This more flexible algorithm has the advantage that it can possibly 
benefit from smaller placing triangulations. The; best strategy for most cases 
seems to be to permute the input randomly rather than to sort it. 

We omit a discussion of the data structures required in the implementa- 
tion. Instead we want to to spend a few words on the arithmetic to be used. 
In a principal way, convex hull computations makes sense over any (compu- 
tationally feasible) ordered field. Natural choices certainly include the field 
of rational numbers as well as certain (or even all) algebraic extensions."'^ 
In the following we focus on the rationals and their extensions by radical 
expressions. 

Examining the Algorithms A and B more closely reveals that computa- 
tions within the coordinate domain are necessary only to decide whether a 
given point is contained in a given afSne hyperplane and, if not, on which side 
it lies. That is to say, if we are only interested in combinatorially correct out- 
put, it suffices to evaluate signs of determinants without ever knowing their 
precise values. This approach, sometimes called robust or exact geometric 
computation, is taken in several computational geometry libraries including 
CGAL [10] and LEDA [23]. It has the advantage that it can be extended to 

^ For a more general perspective on the subject see the book of Blum et al. [5]. 
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radical extensions of Q rather easily. This is implemented in LEDA and in the 
CORE Hbrary [29]. 

If, on the other hand, we do need the exact facet normal vectors then 
we have to use an exact implementation of the arithmetic functions of our 
coordinate domain. Today's standard for long integer and rational arithmetic 
is the GNU Multiprecision Librajry (GMP) [17]. In fact, this is the arith- 
metic used in polymake's implementation of Beneath-and-Beyond as well as 
in the codes of Avis [2] and Fukuda [11]. There are packages in computer 
algebra and computational number theory which can perform computations 
in (arbitrary) finite extensions of Q, but currently this functionality does not 
seem to be available as a stand-alone library which encapsulates a given field 
as a number type to be used in standard C or C-h- 1- code. 

Even if one agrees to compute with rational coordinates, on a technical 
level there is still one choice to be made: One can cither compute with repre- 
sentations of rational numbers directly, or one can translate everything into 
integer coordinates, essentially by scaling. While this usually does not make 
much of a difference for the algorithms, the actual numbers that occur dur- 
ing the computations are different and, in particular, of different sizes. This 
may affect the performance of an algorithm, but often one does not know in 
advance which method is superior to the other. 

It would be interesting to combine the two techniques, exact geometric 
computation and exact coordinates, for the Beneath-and-Beyond algorithm 
in the following way: Use exact geometric computation to produce the com- 
binatorially triangulation, and only in the end compute the facet normals by 
solving systems of linear equations defined by the (d— l)-faces on the bound- 
ary. Since in many cases the arithmetic consumes most of the running time 
it should be possible to save some time this way. The author is not aware of 
an implementation of any convex hull code taking this approach. 

5 Empirical Results 

The analysis in Section 3 showed that, from a theoretical point of view, 
Beneath-and-Beyond seems to be a weak algorithm: It has a super-polynomial 
total running time even on simplicial polytopcs. In this section wc display 
a few computational results, where the performance of polymake ' s imple- 
mentation beneath_beyond is compared to two other programs: Fukuda's 
odd, which implements (dual) Fouricr-Motzkin elimination, sec Zicglcr [30, 
Sections 1.2 and 1.3] or Fukuda and Prodon [13], and Avis' Irs, which uses 
reverse search, see Avis and Fukuda [4,1]. 

The computational results below arc intended to complement the corre- 
sponding data in the paper by Avis, Bremner, and Seidel [3]. The performance 
comparison among the three programs is fair in the sense that all programs 
use the same; implementation of exact rational arithmetic, namc^Iy from the 
GNU Multiprecision Library (GMP) [17]. However, beneath_beyond and 
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cdd both use rational coordinates while Irs uses integers only. This could con- 
tribute to beneath_beyond's superiority over Irs for the "random spheres" , 
see Figures 6 and 7; but this needs further investigation. 

Roughly speaking the input pol3ftopes come in two groups: The first group 
consists of the intricate polytopes discussed in the previous section. These are 

the dwarfed cubes, products of simpliccs, and dwarfed products of polygons. 
While they are known to be computationally hard for all iterative convex hull 
algorithms, this study shows that beneath_beyond performs particularly bad. 
The second class of polytopes investigated arc convex hulls of (uniformly 
distributed) random points on the unit sphere, sometimes called "random 
spheres". Such polytopes are almost always simplicial (in fact, all the tested 
ones were simplicial). The performance of beneath_beyond is clearly better 
than that of Irs, which in turn clearly beats cdd. 

The "random spheres" were constructed with polymake's client program 
rand_sphere which produces uniformly distributed random points on the 
unit sphere with double (i.e., IEEE 64-Bit floating point) coordinates rounded 
to six digits after the decimal point, which were then transformed into exact 
rational numbers; the numerators and denominators of such rational numbers 
typically have 15 to 20 decimal digits each. 

While in [3] the authors also compared various insertion strategies for the 
incremental algorithms we did not do that here: cdd was run with the lexmin 
insertion rule (which is its default behavior). This seems justified in view of 
the computational results in loc. cit., which showed that lexmin is optimal 
among the strategies tested in most cases. 

All timings are taken on several identical Linux machines with an Athlon 
XP 1800+ processor (1533.433 MHz clock, 3060.53 bogomips) and 512 MB 
main memory running RedHat Linux 7.3. The timings were taken via Perl's 
timeO function, and we counted the time in user mode only. The parame- 
ter datasize was limited to 400 MB. Since all programs were run through 
polymake, by calling the clients cdd_ch_client and lrs_ch_client, respec- 
tively, there is a certain additional overhead (due to socket communication 
and data conversion) which should approximately be the same for the three 
programs. This makes up for the slight difi^erence in shape of some of the 
curves as compared to the corresponding ones in [3]. In particular, timings 
below one second are almost impossible to interpret this way. Moreover, the 
implementations of cdd and Irs tested here are more advanced than the ones 
tested in loc. cit. In particular, the old versions of cdd and Irs did not use 
the (considerably faster) GMP arithmetic. 

Due to the fact that the timings include some input/output overhead it 
is impossible to reproduce exact timing data. Therefore we took the average 
over several runs. Additionally, we randomly permuted the order of the input. 
While cdd and Irs arc almost insensitive to the ordering, this can make a 
considerable difference for beneath_beyond. The highest variation among the 
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individual timings of beneath_beyond occurred for the dwarfed products of 
polygons: A more detailed statistical analysis is given in Table 1. 
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Fig. 3. Dwarfed products (;io(s) of five s-gons. In this test s is always oven. Com- 
pare [3, Figure 5]. cdd and Irs: average over 10 runs, beneath_beyond: average over 
50 runs. Memory overflow (more than 400 MB required) in beneath_beyond for 
s > 12. For more details of the beneath_beyond timings see Table 1. 



Table 1. The timings (in seconds) for the individual runs of beneath_beyond for 
the dwarfed products of polygons flio(s), see Figure 3, vary rather strongly. 50 runs 
with a random insertion order were performed. 



s 


average 


minimum 


maximum 


standard deviation 


4 


0.650 


0.35 


1.05 


0.1628 


6 


5.530 


0.78 


23.20 


4.6607 


8 


44.091 


3.26 


226.55 


40.7324 


10 


476.051 


10.00 


1984.59 


452.4177 


12 


1883.964 


239.46 


8560.20 


1803.4397 
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Fig. 4. Dwarfed d-cubes cj, for d even. Compare [3, Figure 4]. cdd and Irs: average 
over 10 runs, beneath_beyond: average over 50 runs, memory overflow (more than 
400 MB required) for d > 16. 



For the "random spheres" we changed our experimental set-up shghtly. 
For each set of parameters d € {3, 4, 5, 6} and s € {100, 120, 140, . . . , 500} 
we produced 10 polytopes as the convex hulls of s random points on the 
unit sphere S'*"^ C R'*. Each of the three convex hull codes was run once on 
each polytope (for beneath_beyond we inserted the points in the same order 
as they were produced). The charts show the average values for each code 
taken over the 10 samples. The individual timings varied only by a little: 
The highest deviations occurring for cdd and d = 6, see the Table 2 for more 
details. 

6 Concluding Remarks 

Authors in discrete geometry often talk about "triangulations of finite point 
sets" . They refer to triangulations of the convex hull of these points. For more 
details see Pfeifle and Rambau [26]. 

Finite simplicial complexes, as combinatorial abstractions of triangula- 
tions, are an indispensible tool in topology. At first sight it may seem acciden- 
tal that notions from topology appear in questions concerning the complexity 
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Fig. 5. Products of two d-simplices. Compare [3, Figure 2]. Average over 10 runs 
for eacii program. berLeath_beyond: Memory overflow (more than 400 MB required) 
for d > 11. 



Table 2. The timings (in seconds) for the individual runs of cdd for s random 
points on the unit sphere in R®, see Figure 7, vary in a rather narrow range. We 
checked 10 cases for each value of s, produced independently at random. 



s 


average 


minimum 


maximum 


standard deviation 


100 


258.997 


247.05 


273.58 


7.5541 


120 


427.075 


416.67 


446.39 


8.3978 


140 


641.947 


627.80 


666.61 


12.0566 


160 


907.759 


895.52 


938.18 


12.9157 


180 


1238.253 


1202.69 


1265.81 


18.9353 


200 


1630.816 


1591.89 


1670.96 


23.6467 


220 


2078.458 


2028.65 


2132.95 


29.9611 


240 


2582.256 


2495.02 


2650.31 


41.1209 




Fig. 6. "Random spheres" with n vertices in dimensions 3 (top) and 4 (bottom). 
Average over 10 polytopcs, each program run only once. Note that our timings 
below 1 second are not very accurate. 
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Fig. 7. "Random spheres" with n vertices in dimensions 5 (top) and 6 (bottom). 
Average over 10 polytopcs, each program nm only oucc. cdd not tested for input 
with more than 240 vertices since it takes about three hours per test. 
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of convex hull computations. However, as developed in the paper [21] there is 
a deeper connection: Up to polynomial equivalence a convex hull computation 
can be replaced by suitable simplicial homology computations. For details on 
algorithms and implementations to compute homology see Dumas et al. [8]. 
The paper [22] by Kaibel and Pfetsch contains more information about the 
complexity status of the convex hull problem. 
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